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Abstract 

We study the evolutionary dynamics of a haploid population of infinite size 
recombining with a probability r in a two locus model. Starting from a low 
fitness locus, the population is evolved under mutation, selection and recom- 
bination until a finite fraction of the population reaches the fittest locus. An 
analytical method is developed to calculate the fixation time T to the fittest 
locus for various choices of epistasis. We find that (1) for negative epistasis, 
T decreases slowly for small r but decays fast at larger r (2) for positive 
epistasis, T increases linearly for small r and mildly for large r (3) for com- 
pensatory mutation, T diverges as a power law with logarithmic corrections 
as the recombination fraction approaches a critical value. Our calculations 
are seen to be in good agreement with the exact numerical results. 
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1. Introduction 



Sexual reproduction is ubiquitous in nature - most eukaryotes reproduce 
sexually and genetic mixing is common in some bacterias 2l|, |25|| . However, 
asexual reproduction is not entirely absent. Microbes such as virus and 



bacteria reproduce asexually most of the time, ancient asexuals which 



have remained exclusively asexual for millions of years persist and human 



mitochondrial DNA has not recombined for a few million years [17|]. It is 
then natural to ask: under what conditions is one or the other mode of 
reproduction preferred? 

A detailed study of theoretical models has been helpful in identifying 
some relevant parameters and conditions. A parameter which plays a cru- 
cial role in the evolution of sex and recombination is epistatic interaction 
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amongst gene loci |16|. Experiments have shown that the individual locus 
do not always contribute independently to the fitness of the whole sequence 



28l . \24i\ and the deviation of the fitness from the independent loci model is 



a measure of the epistatic interactions. The nature of epistasis is important 
in determining whether a mode of reproduction may be viable. For instance, 
in the absence of back mutations and recombination, a finite asexual pop- 
ulation evolving on a nonepistatic fitness landscape accumulates deleterious 
mutations irreversibly (Muller's ratchet) [l^, 0]- But the degeneration can 
be effectively halted if synergistic epistasis is present [l3|, [lOj. On complex 
multipeaked fitness landscapes that incorporate sign epistasis [27|, the effect 
of sex has been seen to depend on the detailed topography of the fitness 
landscape 14, 0|. 

The epistatic interactions play an important role in infinitely large popu- 
lations as well. In a two locus model with the four possible sequences denoted 
by ab, Ab, aB and AB and with respective fitnesses wi, W2, W2, W4(> wi, W2), 
recombination reduces the frequency of the favorable mutant AB when epis- 
tasis parameter e = w^Wi — is positive but increases the AB frequency 
for negative e jsj. In this article, we ask: in an infinitely large recombining 
population if all the population is initially located at the sequence ab, how 
much time T does it take to get fixed to the double mutant AB with fitness 
W4 > Wi7 The fixation time T is expected to decrease with recombination for 
negative epistasis and increase for positive epistasis 20, 6[. These qualitative 
trends are understandable from the results of j^: for e < 0, as recombination 
acts in favor of the double mutant, it will get fixed faster than in the asexual 
case while the reverse holds for e > case. 

The main purpose of this article is to find analytical expressions for the 
fixation time T. To this end, we develop a new method to handle the inher- 
ently nonlinear equations obeyed by the genotype frequencies in the presence 
of recombination (see Section [2]). The basic idea of our approach is that at 
any instant, only one of the genotypes dominate so that the equations can be 
expanded perturbatively in powers of the ratio of the non-dominant genotype 
frequency to the dominant one. 

The rest of the article is organised as follows. We first define the model 
under consideration in Section [2l The dynamics of the population frequencies 
for various choices of epistasis are discussed in Section [31 The fixation time 
defined as the time at which the double mutant frequency reaches a given 
finite fraction is calculated in Section HI The effect of initial conditions on 
fixation time is considered in Section O The last section discusses our results 
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51 which are summarised in Table [H 



52 2. Model 

We consider a two locus model with sequences denoted by ab, Ab, aB and 
AB and respective fitnesses wi,W2,W'i,Wi{> wi,W2,W3). The population at 
these sequences evolves according to mutation, selection and recombination 
dynamics. In such models, several schemes have been used to implement 
these basic processes such as recombination followed by mutation and then 



selection [6|, selection, mutation and then recombination [15[ and selection, 
mutation and recombination appearing as additive terms in continuous time 
models Here we work with a discrete time, two locus model in which 



the mutation occurs after recombination and selection [26|. The mutation 
probability from a to A and 6 to S is given by /i but the back mutations are 
neglected. The recombination between ab and AB or aB and Ab occurs with 
a probability r. Denoting the population fraction at sequences ab, Ab, aB 
and AB at time t by xi, ...,X4 respectively and time t + 1 by x[, ...,x\, the 
time evolution occurs according to the following nonlinear coupled equations: 

/ _ (1 -rP) 
~ 

, _ /i(l - /x)(wiXi - rP) + (1 - fi){w2X2 + rP) 

" wit) 
, _ - \i){wxXx - rP) + (1 - fi){w3X3 + rP) 

, _ fi^jwiXi - rP) + fi{w2X2 + rP) + fijwsXs + rP) + (W4X4 - rP\ ^^ 

w{t) 

Here P(t) = {wiW4Xi(t)x4(t) — W2WsX2it)x3{t))/w^{t) is the linkage dise- 
quilibrium at time t and w{t) = Ylt=i''^kXk{t) is the average fitness of the 
population. 

In the following, we will work with Wi = 1,W2 = W3 and > Wi,W2 
and initial condition Xfc(O) = 6k,i- As a consequence, X2(t) = xs{t) for all 
t > 0. We define the epistasis parameter e = W1W4 — W2 and will discuss four 
separate cases: (i) zero epistasis which requires W2 > 1 as 1^4 > 1 (ii) negative 
epistasis (iii) positive epistasis and W2 > I (iv) positive epistasis and W2 < I 
(compensatory mutation). It is useful to write Xkit) = Zk{t)/ J2^=i ^ji't) 
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where the z^s satisfy the following condition: 



4 4 
i=l 



i=l 



The unnormalised populations z^s obey the following set of equations, 



- n) (^Zi - rD^ + {1 - fi) (w2Z2 + rD^ 



4 = yzi - rDj + 2/i \^W2Z2 + rDj + yw4Z4 - rD 
with the initial condition Zk{0) = 5^1. In the above equations, 



(5) 

(6) 
(7) 

(8) 



D 



Zi 



W4Z1Z4 — w\z\ 



WiZiY 



(9) 



For r = 0, the above model reduces to the standard quasispecies model 



equations 11 


. For 








^ 1 




Z2{t) 


^ /i 




Z4{t) 





■^2 



1(72 + 1 tf. 



2w2 wi — wi 



W2 — 1 W4 — W2 W2 — I W4 — 1 



(10) 

(11) 

(12) 



56 
57 
5B 
59 
60 
61 



Figure [T] shows the time evolution of the populations ^^'s when r = for 
negative and positive epistasis. 

With nonzero recombination, it does not seem possible to solve the above 
equations for Zi{t) exactly due to the bilinear terms in linkage disequilbrium 
D. However, in the next section, we will obtain approximate expressions for 
the unnormalised population Zi. 



62 3. Time evolution of populations 

63 As we shall see, the dynamics of population Zj's can be divided in following 

64 three dynamical phases: (i) zi ^ Z2,Z4 (phase 1) (ii) Z2 ^ 2^1,-24 (phase 
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Figure 1: Recombination probability r = 0: Time evolution of zi (solid), Z2 (broken) and 
Z4 (dotted) for ^2 = 2, W4 = 3, /i = 10~^, and zi (+), Z2 (x) and Z4 (□) for t«2 =2,W4 = 
5,fi= 10^^ using exact equations ©-(IS]). 
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II) and (iii) 2:4 ^ -21,-22 (phase III). Thus we can expand equations ([H])-(IHD 
in powers of Z2I Zx^z^j Z\ in phase I, zi/ Z2, z^j Z2 in phase II and similarly, 
Zi / Z2 / za in phase III. The time scale at which a phase ends is obtained 
by matching the solutions of the relevant populations in the two phases. 
In the crossover region however the above assumptions are not expected 
to hold strictly. But as we shall see, the fixation time is nevertheless well 
approximated. Note that the perturbation expansions here are different from 



a small r expansion [15 
3.1. No epistasis 

When W4 = w^, the epistasis parameter e = 0. The dynamics of popula- 
tions Xj's and -Zj's evolving for such a fitness choice is shown in Fig. [2J For 
e = 0, the linkage disequilibrium D{t) obeys the following evolution equation: 



D' 



rDSiS2] 



oc 



[(1-^- 

^2 (-21-^4 



fiW2yzi + 2^2(1 - fi + fiW2)z2 + W2Z4 - r(l - /i)^(l - ^2)^-05*1] 



- rDSiS2 



(13) 



where Si = J2t=iZi and S2 = J2t=i'^iZi- As D oc ^12:4 — for e = 
, it follows that D' ~ D. Therefore, if the population has zero linkage 
disequilibrium to start with, it remains zero for all times in the absence 
of epistasis Thus for e = 0, the population ZkS obey following linear 
equations, 

z[ = (1 - fifzi 

Z'2 = jU(l - fl)Zi + (1 - fx)w2Z2 
z'^ = ^^Zi + 2fiW2Z2 + WaZa 

which can be easily solved. For small /i, we obtain 



Zi{t) 
Z2{t) 

Z^{t) 



1 



1 



W2- 1 











\W2 





(14) 
(15) 

(16) 



75 From the above solution, we see that at time ti at which phase I ends, 

76 ziZa — Z2 = Z4 — Z2 = so that Z4 = Z2 and therefore the phase II is absent 

77 in this case. 
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Figure 2: Non-cpistatic interaction: Time evolution of z\ (solid), Z2 (broken) and 04 
(dotted) for wi = 1.25, i(;4 = 1.5625,7' — 0.1, /i — 10^^ using exact equations ([H])-®. The 
normalised fractions are shown in the inset. 
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Figure 3: Negative epistasis: Time evolution of zi (solid), Z2 (broken) and Z4 (dotted) for 
W2 = 2,W4 — 3,7- = 0.1, /I = 10^^ using exact equations ([6])-([8]). 



3.2. Negative epistasis 

We now consider the case when epistasis is negative, < W2. The time 
evolution of populations for this fitness scheme is shown in Fig. [1] for r = 
and Fig. [3] for r > 0. In both cases, the population zi dominates at short 
times followed by Z2 and finally Z4 takes over. This behavior is also reflected 
in the normalised populations x^s shown in the inset of Fig. [3l 
Phase I. Since ^^(0) = 5^,1, initially Z2,Z4 ^ zi so that the sum J2t=i ~ -^i- 
Using this in the expression for D, we find that D ^ ^42:4 — Taking 
/i — >• in ([6])- ([8]), we then obtain 



z'l ~ Zi + r 



f W2Z2 — W4Z4Z1 \ 
\ J 

z'2 ^ /i^i + ^2-22 - r 



w\z\ — W4Z4Z1 



Zl 

^ /i^Zi + 2^W2Z2 + W4Z4 + r 



W2Z2 — W4Z4ZI \ 



Zl J 

Since Z2/ zi^ -24/^1 <^ 1, we can write z[ ~ Zi which immediately gives zi{t) 
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1 for t < Ti where Ti is the time at which phase I ends. In the equation for Z2 
(2:4) , the first term on the RHS is the mutation term which signifies that even 
if W2 {W4) and r are equal to zero, the population Z2 {z^j will remain nonzero 
due to a constant mutational supply from ab. The second term on the RHS 
of the Z2 equation is the selection term and the last two terms are due to 
recombination. It turns out that the recombination term can be ignored to 
yield Z2 ~ fiZi + W2Z2- This approximation is later justified by showing that 
indeed the last two terms are negligible compared to W2Z2- In summary, we 
have 



fIZi + W2Z2 



.2^2 



~ fi Zl + 2fiW2Z2 + (1 — rjw^Zi + rw2 — 
On solving the above equations, we obtain 



Zl 



(17) 
(18) 

(19) 



Zi{t) 


^ 1 


^2(t) 


^ /i 


z^it) 





Wr, 



1 



W2 - 1, 

(■^4(1 — r))* 



1 



+ 



2fi^W2 



+ 



1 1 

rfi W2 



W4[l — r) — 1 

(^4(1 — r))* — W2 



W2- I 

t „,,2t 



{wi{l — r)y — w\ {wa{'^ — r)y 
^4(1 - 
.(^4(1 



(20) 
(21) 
- 1 



r) - W2 
-r)y 



^2 , (w4(i - r)y -4; ^ 



{w2 — ly I ^4(1 — r)— ^2 W4{l — r)—W2 W4 

We check that the r = limit is recovered from the above solution. From 
the last equation, we find that for r > 0, the growth rate of population 24 is 
given by max{w4(l — r), Since r must be positive, for e < 0, the growth 

'4 



rate of z^ is Wn and we have 



Z4{t) 



1 1 

rw2 



{rw4 + \e\){w2 - ly 



X wf 



r > 



(23) 



84 Using the above solutions, it can be checked that W2Z2 ^ r(w|z| — W4Z4) 

85 thus justifying (fT8l). 

It is evident from (12T1) and (1231) that when Z2 becomes one, Z4 = rwl/ (rw4+ 
|e|) < 1 so that Z2 intersects Zi before Z4. The time ti at which 22(^1) = 1 is 
given by 

r,^^Xn{^) (24) 



In 



/i 
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86 which is independent of r. 

Phase II. After time ri, the population Z2 ^ -21,-24 and therefore the sum 
J2i=i^i ~ 2^2- For weak selection, this gives D ~ —Z2/2. 



phase, ZkS obey the following equations: 



4 



^2 



4 



r 

^1 + 



Thus in this 

(25) 
(26) 
(27) 



As the equation for Z2 is decoupled from zi and Z4, we can first solve for Z2 
and then use the solution to find zi and Z4. This finally gives 

f \ t — Tl 

2, 



Z2{t) 
ZA{t) 
Z,(t) 



^2 Ti, 



,t-Tl 



t-T, / N r[w2- {r/2)] ''^ -W4^ , , 
2 W2- (r/2) - W4 



. ^ r [W2 - (r/2)]*-- - 
'^^"^^ + 2 «;2-(r/2)-l 



-^2(ri) 



(28) 
(29) 

(30) 



87 where -2i(ti) = Z2{ti) = 1, ^4(ri) = rw|/(ru»4 + |e|). 
The time T2 at which 2^4 overtakes Z2 is given by 



T2 = Ti + 



ln[2^/;4/(2«;2-r)] 



In 



2(rw4 + |e|)(w2 - W4 



[2(w2 — W4) — r] rw 



re 



r > 
(31) 

We check that for e — 0, the time T2 — ti during which phase II is present 
vanishes. To understand how T2 varies with r, consider the r-dependent term 
g{r) in T2 which can be rewritten as 



In [2^4/(2^2 - r)] 



In 



1 + 



ln(tt;4/w2) 
1 

In(w4/w2) 



In 



1 + 



e| (1 - r)(r + 2(w4 - W2)) 
r |e| + rw4 — W2 [2(^2 — ^4) — r] 
(1 -r)(r + 2(^4-^2)) 



r |e| + ri(;4 — W2 [2(^2 — ^4) — r] 



In [1 + 7^] 



The ratio 7^ = 1 when r satisfies the quadratic equation Wgr^ + (?i;4 — 
i(;2)(2w| — 'u;4)r — {wi — W2)\e\ — {r — r+)(r + |r_|) = where r+(r_) is 
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the positive (negative) root of the quadratic equation. For r ^ r+, the ratio 
7^ > 1 so that ln(l + 7^) ^ ln7^. For r > r+, 7^ < 1 and ln(l + 7^) ^ 7^. 
Using these approximations in the expression for g{r) above, we find that 

^(r)~P^[^^'^^~^^^/"] (32) 
1 |e|/r , r > r+ 

Thus the time T2 — ti decreases slowly as ln(l/r) for small r and as 1/r 
for large r. Due to these properties of T2, the single mutant population can 
dominate for appreciable time interval for small r. 

Phase III. For t > T2, the population z^':^ Zi, Z2 so that J2i Zi ^ Z4,. Due to 
this, D ~ (zi/wi) — {{w2zl) / {wizi)) . For small fi, the equation for z^s in 
(El)- ([8]) can thus be simplified to give 

z'^ ^ W4Z4 (33) 
4 ^ W2Z2 (34) 
, , r \ rw^z^ 



z[ ^ 1 \z, + ^^ (35) 

V W4J wi Zi 

where we have neglected the recombination term contribution to the equation 
for Z2 by assuming ^22:2 ^ rzi/w^ (see below). From the first two equations, 
we see that Z2 ~ w\ and 2:4 ~ W4. Thus the last term in the equation for zi 
can contribute when e < 0. Explicitly, we obtain 

z^{t) = w'r'zi{T2) (36) 

Z2{t) = w\~^'Z2{t2) (37) 

zAt) - (1 - -r-.ir2) + -f M)/-;; - 

Wa wi Zi{T2) (1 - (r/w4)) - {W^/Wi) 



91 where Zk{T2) are given by (l28l) -(|30ll at time t = T2. From the above solution, 

92 it is easily verified that W2Z2 ^ rzi/w4 is a good approximation for t > T2. 

93 5. 5. Positive epistasis 

94 We now turn to the case when epistasis is positive. The condition W4 > w| 

95 can be satisfied for W2 < I and W2 > I. For W2 > I, the time evolution of 

96 populations for this fitness scheme is shown in Fig.|l]for e < rw^ and e > rw^. 



11 



97 The reason for this distinction will be explained below. The dynamics of the 
9B populations Zi, Zi and 2:4 for W2 < 1 are shown in Fig. O Note that phase II 

99 is absent in all these cases. 

Phase I. As discussed for e < 0, in this phase, zi ^ Z2, Z4 and can be well 
approximated by one, Zi{t) ~ 1 for t < ri. Then the populations Z2 and 2:4 
obey the following equations: 

Z2 ~ fi + W2Z2 + r{w4Z4-wlz2) (40) 

Z^ ^ fl + 2fIW2Z2 + W4{1 — r)Z4 + rW2Z2 (41) 

Figures m and [5] show that initially Z2 > z^ but after some time (< ri), Z4 can 
overtake Z2 while both Z2,Z4 < 1. This behavior is characteristic of positive 
epistasis as can be seen from Fig. [T] for r = also. For this reason, the 
population Z2 can get a contribution from Z4 in phase I when e > and we 
need to retain the r-dependent term in the equation for Z2. The equation 
for Z4 remains the same as for negative epistasis. Since Z4 appears linearly 
in the above equations for Z2 and Z4, it is possible to eliminate Z4 from the 
equation for Z2 and express it in terms of Z2 alone. This gives a three term 
recursion relation for Z2: 

Z2{t +1) = (1 + rw4, — W4)[i + rw4[j? + {2iJ,rw4 — 104, + rw4)w2Z2{t — 1) 
+ {w2 + W4 — rw4)z2{t) + rw4wlzl(t — 1) — rwlzl{t) , t > ^42) 

100 with initial conditions -22(0) = and -22(1) = f^- We will find the solution 

101 to the nonlinear equation for Z2{t) iteratively [1]. We first find the solution 

102 fo{t) of the above difference equation for Z2 when the nonlinear terms are 

103 set to zero. The corrections to Z2{t) arising due to nonlinearity will then be 

104 determined by writing Z2(t) = /o(t)(l + fi(t)). 

The solution fo{t) satisfies the following linear, inhomogeneous difference 
equation: 

/o(t + 1) = So/o(t) + Co/o(t - 1) + Ao , t > 1 (43) 

where Aq = {l+rw4-~W4)fi+rw4fi'^, Bq = {w2+W4 — rw4),Co = {2firw4 — W4 + 
rw4)w2- The solution of this linear equation subject to /o(0) = 0, /i(l) = A* 
can be found by using the method of variation of parameters [ij and is given 
by 

= /^(l - «+) - ^ /i(l - a-) - Aq ^ Aq 

(a_|_ — «_)(1 — «+) ^ (a_ — a+)(l — a_) ^ (1 — «_)(1 — a+) 
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Figure 4: Positive epistasis: Time evolution of zi (solid), Z2 (broken) and Z4 (dotted) for (a) 
W2 = 1.25,^4 = 1.8125, r = 0.1, = 10"*^ (b) W2 = 1.25, W4 = 1.8125, r = 0.4, = 10"^ 
using exact equations ([51)- (HI)- The normalised fractions are shown in the inset. 



where 



{w4 — rw4 + W2) ± a/ (^4 — rwi — W2Y + 8j2rw4W2 



105 
106 
107 
108 
109 
110 
111 



For our purposes, it is sufficient to retain terms to in the last expression 

which gives 

/i^rw4 [2^2(1 - W2) + (1 + W2){w4 -W2- rw^)] 



{W2 - 1)(W4 



+ 



fi'^rw 4(104 + W2 — rw4) 
(w4 — rw4 — W2y 



[W4 — rw4 



W2 — rw4 
*- 1" 



Wo 



W2 



W4 — rw4 — 1 



- 1 

(44) 



Note that the above solution consists of two growth rates for Z2 namely W2 
and 104(1 — r). 

Using Z2{t) ^ /o(^) in ( 14TI) and keeping terms to we get ( !22l) for 

Z4(t). It follows that the population Z4 does not grow if r > Tc = {W4 — l)/w4 
for W2 < 1- But for 1^2 > 1, the population Z4 always grows and the growth 
rate is given by max{t/;4(l — r),w\} as in the case of negative e. In the 
following, we will discuss the two cases W2> 1 and 1^2 < 1 separately. 
W2 > 1: When W2 > I, due to ( !22l) . the following subcases arise for Z4{t): 



Z4{t) 



/i^e(l— r)(iii4— ri04+ui2) 
(e—rim){w4—rwi—l){w4,—rwi—W2) 

{rW4—e){w2 — l)''^ 2 



X (w4 — rw4Y ) ''^ < e/w4 
, r > e/w4 



(45) 
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For r > e/w4, it follows from (HSj) that when Z2 becomes one, = rwl/iTW^^— 
e) > 1 so that ^4 hits unity before Z2 and thus phase II is absent. The time 
Ti at which ,24 (ri) = 1 is given by 



In W2 



In 



(— — ^ ) + 777- — In (^4 - -) , > e/wi 



(46) 



112 For r ^ e, the last term in the above expression (and hence ti) increases as 

113 ~ —e/{rw/i) with increasing r. 

For r < e/wi, using fH5|) . we find that the time ri at which Z4(ti) = 1 is 
given by 



ri(r) 



1 



ln(w4 — rwi 
For r = 0, we have 



In 



(e — rw4)(ti?4 — 1 — rwi){w4^ — W2 — rw4) 
jj?e{l — r){wi + W2 — rw4) 



n(0) 



In W4 



In 



(W4 - 1)(W4 - W^2) 
/i2(w4 + W2) 



r < e/w4 
(47) 

(48) 



which matches the one obtained using (fT2|) or (!22|) for W2 -C W4. To find the 
behavior of ri for r -C e/w^, we rewrite the expression for Ti(r) as 



ln(ti)4 — rw4)Ti{r) — lnw4ri(0) 



In 



e( 1 — r J 



rw4 
W4 — 1 



1 - 



rw4 



W4 — W2 



rw4 



Wi + W2 



Using the inequality r < rw^ < e < W4 — W2 < W4 — 1 < + W2 iia the last 
equation, we find 



ln(w4 — rw4)Ti{r) ^ lnw4ri(0) + In I 1 



rwn 



The above expression can be further simplified to give 



ri(r) ^ri(O) 1 + 



In W4 



(49) 



(50) 



114 which shows that Ti increases linearly with r for r <^ e/w^^. 

115 W2 <!'■ It is known that for an infinite population, there exits a critical 

116 recombination fraction Tc beyond which a population initially located at ah 
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Figure 5: Compensatory mutation: Time evolution of zi (solid), Z2 (broken) and Z4, 



(dotted) for ^2 = 0.75, W4 = 1.25, r = 0.175, ^l ■ 
normalised fractions are shown in the inset. 



10 using exact equations 



The 
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cannot cross the intermediate fitness valley and reach the double mutant 
fitness peak For our model, as discussed above, Z4 can grow 

(and hence can be fixed) provided ^4(1 — r) > 1 or r < rc ~ (^4 — l)/w4. 

Due to fH4|) obtained by dropping nonlinear terms in the equation for Z2, 
both Z2 and 2:4 increase as {w^^ — rw4^y. But Fig. [5] shows that Z2 and Z4 do 
not continue to grow in this manner but rise sharply as the end of phase 
I approaches. To understand this, it is essential to include the nonlinear 
terms in the equations (HOj) and fHTl) for Z2 and Z4. To this end, we write 
Z2{t) = fo{t) [1 + fi{t)] where /o(t) given by dH]) reduces to 



, n'^rw4{w4 + W2-rWi) 
1 — W2 {W4 — rw4 — W2) 



(w4 — rw^y 
Wi — rw4 — 



(51) 



for W2 < 1,^4(1 — r) > 1. As we shall later, fi{t) remains close to zero for 
short times but contributes substantially at large times. Using this form of Z2 
in (H2I) and neglecting the quadratic terms in /i, we find that fi{t) obeys the 
following approximate linear inhomogeneous equation with time-dependent 
coefficients, 

hit + 1) = B^{t)f^{t) + Ci{t)h{t - 1) + A^it) (52) 

where the coefficient 

Bi{t) = W2 + — rw4 — 2rwlfo{t) (53) 
Ci{t) = {2firw2 — + rw4)w2 + 2rw4wlfo(t) (54) 
A(t) = r{w4-l)wlfo{t) (55) 

For — rwi > 1, we define e = Tc — r. For e — > 0, fo(t) ~ a(l + 6^4)* where 

/iVc(l + W2) 



[1-W2ye 



(56) 



and the coefficients are given by 



B,{t) ^ W2 + l-2r,wlfo{t) (57) 
Ci(t) ^ ~W2 + 2nw4wlfoit) (58) 
A,{t) ^ rlw^wlfoit) (59) 



16 



Writing the difference equation (132]) for as a differential equation, 
we get 

dh it) 1 - B^jt) - Crjt) 

with the initial condition /i(0) = 0. It is straightforward to solve the above 
differential equation and we obtain 

/iW«5[('' + ''(l+«"')T-l] (61) 



= , ..:,z. m 



where 

1 -W2 

1 — W2 + 2rcWiW2a 
c = 1-b (63) 

a = — (64) 

The second term in the parentheses in the above equation can be neglected 
for t <^ ti = In [(1 — i(;2)/(2rcW4w|a)] /ew4 and since a ~ fi^, we obtain 
/i(t) ~ below this time scale. For larger times t ^ ti, the first term in 
the parentheses can be ignored and for e — 0, we obtain /i(t) ~ e''''*. Thus 
Z2(t) = fo{t){l + /i(t)) increases as 

120 Thus at times close to the end of phase I, Z2 increases at a faster rate. 

To find the time ti at which phase I ends, we first calculate 2:4 (t) using 
Z2it) = /o(t)(l + hit)) in m- This yields 

^ /o(^ + !)[! + flit + 1)] - - ^2/o(t) [1 + /i(t)] + rwlf^jt) [1 + 2/1 (t)] 

rw4 

(66) 

On expanding /o(t + 1) and fi(t + 1) for small e, we obtain 

/o(t + l) ^ (l+e«;4)/o(t) (67) 
/i(t + l) ^ /i(t) + iceaw4(l + eii;4)* [6 + c(l + eM;4)*]°"' (68) 
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Substituting this in the above expression for 2:4 (t) and using z^^ti) = 1, we 
find that Ti is determined from the following equation: 



'1 — ^2)0 , ^ , arc(c + 2aw|) 2 



-1/1(1 + y2)H y^y2 = rW4 + H (69) 

where yi = (1 + 6^4)^^ and 1/2 = (6 + c(l + 6^4)"^^)". As this equation is 
difficult to analyse, we consider only the fastest growing term to obtain the 
following approximate equation for ti. 

yly2Ci'i^c [c + 2awl] ~ 2rw4 (70) 

Phase III. For t > ti, the populations ^^'s obey the equations as for 

e < and the corresponding solutions are given by fl36l) - (l38l) with T2 replaced 
by Ti. For W2 > I, Z2 and 2:4 grow exponentially fast with their respective 
fitnesses but zi decays with time. The rate of decline is determined by the 
ratio tf|/(w4 — r). If this ratio is larger than unity, Zi ~ (iy|/w4)* and as 
((^4 — r)/wiY otherwise. For 1^2 < 1, ^2 decays with time while ^4 continues 
to grow. The remarks above for zi behavior when ^2 > 1 hold for 'W2 < ^ 
case also. 



129 4. Fixation time 

As seen in the last section, the unnormalised populations z^'s vary ex- 
ponentially (or faster) with time so that the normalised population X4 will 
reach unity asymptotically. Therefore, we define the fixation time T as the 
time when the population fraction 0:4 (T) = 1 — 5 where (5 — >■ 0. In terms of 
ZkS, this condition gives 

Zi(T) + 2z2{T) - -r^MT) = (71) 

1 — 

130 where Zk{T) in the above equation is the population fraction in the Phase 

131 III at t = T. Another reason why 5 > is that for 5 = 0, the above equation 

132 cannot be satisfied as both zi and Z2 are always positive. 

For r = and 1^2 > 1, since 2:1 ~ 1, we can write dz/^iT) ^ 2z2{T) which 
gives 

2{l-5){w,-W2){w,-l) A 

^5{Wi + W2){w2-l) J 

133 which decreases monotonically as W4 increases. 



T ^ ^ In 

\n{wi/w2) 
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134 ^0 epistasis 

Since z^it) = zlif) due to f|T5l) and (fT6l) . the condition flTTl) simplifies to 



give 2;2(^) = (1 + -2^2 (^))v]—^ which leads to 

(W2 - i)v/r^ 

/x(l-v/r^) 



m W2 



1 + 



(73) 



135 For W2 = 2,^4 = 4, /i = 10~^ and 5 = 0.01, the above expression yields 

136 T = 27.56 in excellent agreement with the result of our exact numerical 

137 iteration which gives the fixation time equal to 28 for various values of r. 

138 4- 2- Negative epistasis 

Using (1361)- (I38D at t = T in (JTlD, we have 

2 / \ T—T2 r / \ T—T2 

™^ ^"^^ 2--^^] =0 (74) 



W4{r — e) 1 — 5 \'"^2 

Since the first term on the LHS is exponentially decaying, we neglect it to 
obtain 

^2(1-5)^ 



6 



(75) 



139 where T2 is given by (131]) . As discussed in Sec. 13. 2[ since T2 decreases with r, 

140 the fixation time T is a decreasing function of recombination probability r 

141 when epistasis is negative. A comparison of the analytical estimate with the 

142 exact numerical result for two sets of parameters shows a good agreement 

143 (see Fig. [6]). 

144 4-3- Positive epistasis 

Dividing both sides of flTTl) by W2 and using flS^ - flSS]) . we find that the 
contribution due to Zi term can be neglected as it is exponentially decaying. 
This gives 

r / \ T—T\ 
I W4 



and hence 



2.2(r.)-^^^) (76) 

^ In [2(1 -5) A] ^ ln^2(ri) 
In(w4/M;2) \n{w4,/w2) 



145 We first consider the W2 > 1 case followed by W2 < 1- 
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Figure 6: Negative epistasis: Fixation time as a function of r obtained using exact iteration 
(•) and analytical results (o) given by (|75|) for two sets of fitnesses and /i — 10~^, S = 0.01. 
The data for W2 — 2.0, W4 = 3.0 have been shifted by a constant by adding 12. 
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Figure 7: Positive epistasis: Fixation time as a function of r obtained using exact iteration 
(•) and analytical result (o and x) given by (|77p and ri (upto a constant) respectively for 
two sets of fitnesses and fi = 10^^, S = 0.01. 



146 W2 > I- In the above expression, ^2(^1) is given by (jUj) and ri by (j47j) for 

147 r < e/w^ and fH6l) for r > e/t(74. Computing T using these formulae in fl77|l . 

148 we obtain the fixation time as a function of r shown in Fig. [7] (open circles) 

149 for fitness choices W2 = 1.25,^4 = 1.8125, e/w4 ~ 0.14 and W2 = 1.25,^4 = 

150 2.5, e/w4 = 0.375. The analytical data are seen to be in good agreement with 

151 the exact numerical results except in the vicinity of r = e/w4. Figure [7] also 

152 shows Ti which displays a similar behavior as T. We have already seen that 

153 the time ti increases linearly with r for r <^ e/w^ but weakly for r ^ e/w^. 

154 Thus the fixation time T for e > 0, 1^2 > 1 increases fast for small r but is 

155 weakly dependent on r for r > e/w^. 

W2 < The inset of Fig. [8] shows that the fixation time diverges as r ap- 
proaches critical recombination probability = {w^ — l)/w4. From fl77|) . 
the fixation time T for W2 < I can be calculated using ti from fl69|) and 
-22(^1) = /o(ti)(1 + /i(ti)). The time T thus obtained (open circles) when 
plotted against Tc — r is compared with the results from exact iteration in 
Fig. [8] and shows a good agreement. The approximate ri obtained using (170!) 
is also plotted which well approximates the fixation time T. Therefore, it is 
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Figure 8: Compensatory mutation: Fixation time as a function of r obtained using exact 
iteration (•) and analytical results (o and x) given by (|77|) and ((7D|) respectively for 
parameters W2 = 0.75, u)4 = 1.25, /i = 10^^,(5 = 0.01 with Tc = (w4 — l)/w4. The solid 
line has a slope equal to —1. 
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sufficient to analyse (170]) in order to understand the behavior of the fixation 
time as r — i> r^. For e ^ 0, fl70l) can be written as 



yh^arc [c + 2awj] ^ 2rcW4 
On taking logarithms both sides, the above equation reduces to 



1 - W2)e 



In- 



^4(1 - W2) 



1 — W2 + rcW4)w2al 



(78) 



(79) 



where we have neglected the linear term in ri as compared to the exponential 
term in ti. Writing a = a/e, we finally obtain 



1 



ew4 



In In 



^4(1 - W2)e^ 



W2 + rcW4)w2a? 



-In 



2r'lw2a 



(l-w;2)eVJ 



(80) 



156 which decays slower than 1/e (see Fig. [8]) due to the logarithmic corrections. 



157 5. Initial Condition with nonzero linkage disequilibrium 

158 So far, we have discussed the population dynamics starting with an initial 

159 condition in which only one genotype has a nonzero population. In this sec- 

160 tion, we consider the situation when a small finite frequency at the interme- 

161 diate loci is also present at t = i.e. xi(0) 7^ 0, 0:2(0) = ^3(0) = (1 — a;i(0))/2. 

162 As the analytical method presented in the last sections assumes that all but 

163 one frequencies is rare at a given time, it seems difficult to obtain analytical 

164 results. Therefore we present numerical results to show how the change in 

165 initial condition affects the fixation time. 

166 As shown in Fig. [9], due to a nonzero population at intermediate loci, the 

167 fixation time at a given r is reduced as compared to the situation when only 

168 the genotype ah is present initially. For negative epistasis (Fig.lH^), the trend 

169 in the generalised situation appears similar to that discussed in Fig. [H] in that 

170 the T decreases slowly for small r but fast for large r. However for positive 

171 epistasis with W2 > 1 (Fig. [9)d), the fixation time remains roughly constant 

172 and unlike Fig.[7]does not increase. For compensatory mutations, the fixation 

173 time shown in Fig. [TO] increases as r approaches a critical recombination rate 

174 and diverges slower than (r^ — r)^^. 
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Figure 9: Fixation time as a function of r obtained using exact iteration (•) when a;i(0) = 
0.95,a::2(0) = a;3(0) = 0.025 for the same parameters as in Figs. [6] and [T] 



175 6. Conclusions 

176 In this article, we have studied the dynamics of a 2 locus model in which 

177 the population evolves deterministically under mutation, selection and re- 

178 combination. As the recombination process makes the equations nonlinear, 

179 in general it is difficult to study such problems analytically. Here we have 

180 developed an analytical method to find the fixation time to the best locus 

181 for various fitness schemes. 

182 The fixation time T is one of the measures for judging whether recombi- 

183 nation is beneficial for a population j^. If the fixation time decreases with 

184 increasing recombination rate r, one may deduce that recombining popu- 

185 lation has an advantage over an asexual one. Our calculations show that 

186 when the epistasis parameter e is negative, the fixation time decreases fast 

187 for large r which suggests that high recombination rate may be beneficial for 

188 populations with negatively epistatic fitness. 

189 In a fitness landscape with positive epistasis and fitness increasing mono- 

190 tonically with mutational distance from the initial sequence, the fixation 

191 time is shown to increase with recombination. As already discussed in the 

192 Introduction, the result that T increases with r for positive epistasis is in 

193 qualitative^ agreement with the expectation from the results of Eshel and 

194 Feldman [[5]. Our analytical calculations show that the functional behavior 

195 of time T depends on the ratio rw^/e. The fixation time is shown to increase 

196 linearly when r -C e/w^ but remains roughly constant for r 3> e/w4. 

197 A compensatory mutation is said to occur when the fitness loss caused by 
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Figure 10: Fixation time as a function of r obtained using exact iteration (•) when xi (0) = 
0.95,2:2(0) = a;3(0) = 0.025 for the same parameters as in Figs. El The numerically 
determined critical recombination rate Tc ~ 0.2002 and the solid hne has a slope equal to 
-1. 
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Title 


Conditions 


Fixation time T 


e = 


r > 


Independent of r (1731) 


e < 0, W2 > 1 


r > e 


Decreases with r (|75ll 


e > 0, W2 > 1 


r < e/w^ 


Increases linearly with r (1501) 




r > e/wi 


Increases weakly with r fj46|) 


e > 0, W2 < 1 


r < Vc 


Diverges with Vc — r flHOl) 




r > Vc 


Infinite 



Table 1: Table summarising the dependence of fixation time T on recombination proba- 
bility r for various choices of epistasis e. 



198 one mutation is remedied by its epistatic interaction with a second mutation 

199 at a different site in the genome. For such a fitness scheme in which the 
initial and final fitness hills are separated by a fitness valley, it is known that 
an infinitely large population cannot cross the intermediate valley beyond 
a critical recombination rate rc js], [^. This implies that the fixation time 
diverges as r approaches Vc- Our exact numerical results for fixation time 
when plotted against — r on a double logarithmic scale indicate a power 

205 law decay. Assuming that the divergence is purely algebraic, a fit to the 
numerical data then gives T ~ (r^ — r)~°'^^. However our calculation that 
takes the nonlinearities into account shows that the divergence is actually a 
power law with logarithmic corrections. 

Here we have focused on the evolution of deterministic population but it 

210 is important to include drift effects as the real populations have a finite size 

211 A^. The finite population problem with compensatory mutation has been 
studied in certain parameter regimes using simulations and within a diffu- 
sion approximation. The analytical calculations of (26l . lo], [sj and numerical 



200 
201 
202 
203 
204 



206 
207 
208 
209 



212 
213 



220 



214 simulations of |18l . |23| for a two-locus model with compensatory mutation 

215 indicates that for fixed s and A^, a finite population manages to reach the 

216 best locus for any r and the fixation time increases with r as seen for in- 

217 finite population. However it is not clear how the population approaches 

218 the infinite A^ limit. For e > 0,W2 > 1, it is found numerically there exits 



219 a critical epistasis value below which the fixation time decreases [20|. An 



analytical understanding of such interesting aspects of the interplay between 

221 recombination and drift remains an open problem. 
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